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Abstract 

We propose a computationally efficient atom-superposition-based method for simulating spin- 
polarized scanning tunneling spectroscopy (SP-STS) on complex magnetic surfaces based on the 
sample and tip electronic structure obtained from first principles. By directly calculating bias- 
integrated terms, we avoid numerical errors of the differentiation of the tunneling current with 
respect to the bias voltage in determining the differential conductance (dl/dV). The capabilities 
of our approach are illustrated for a Cr monolayer on a Ag(lll) surface in a noncollinear mag- 
netic state. We find evidence that the simulated tunneling spectra and magnetic asymmetries 
are sensitive to the tip electronic structure, and we analyze the contributing terms. Related to 
SP-STS experiments, we show how to simulate two-dimensional differential conductance maps 
and qualitatively correct effective spin polarization maps on a constant current contour above a 
magnetic surface. Moreover, we derive alternative expressions for dl/dV and suggest a combined 
experimental-theoretical procedure that might help to extract surface and tip electronic structure 
information from experimentally measured spectra. 

PACS numbers: 72.25.Ba, 68.37.Ef, 71.15.-m, 73.22.-f 
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I. INTRODUCTION 



The scanning tunneling microscope (STM) and its spectroscopic mode (STS) proved to 
be extremely useful for studying local physical and chemical phenomena on surfaces since 
the invention of the STM 30 years ago |2J. The progress of experimental techniques in the 
last two decades was remarkable, thus, more sophisticated theoretical models and simulation 
tools are needed to explain all relevant details of electron tunneling transport measurements 
3j,|4|. STS theory and applications are recently focused on extracting surface local electronic 
properties from experimental differential conductance (dl/dV) data fl-EJ. The role of the 
tip electronic structure has been identified to be crucial on the dl/dV tunneling spectra, 



see e.g. Refs. 



and a theoretical method has been proposed to separate the tip and 



sample contributions to STS 



111- 



An emerging research field in surface science is the investigation of magnetism at the 
nanoscale and atomic scale with the aim of achieving ultrahigh information density for data 

13]. Spin-polarized scanning tunneling microscopy (SP-STM) [l^] 



storage purposes 



12 



is admittedly an important tool for studying magnetism on surfaces. Recent experimen- 
tal advances using this technique allow the investigation of complex magnetic structures 
(frustrated antiferromagnets, spin spirals, skyrmion lattices, etc.) 15Ml8|. Spin-polarized 



scanning tunneling spectroscopy (SP-STS) has recently been used to find inversion of spin 
polarization above magnetic adatoms 19M21]. and the effect has been explained theoretically 
221 ] . Furthermore, SP-STS is useful to study atomic magnetism 23], or many-body effects 



25] as well. The effect of differently magnetized surface 



on substrate-supported adatoms 
regions on SP-STS has also been reported |26|, 12711, a rid the role of tip effects on SP-STS 



28 



29] and on achieving giant magnetic contrast 



30] have also been highlighted 



In the present work, we extend our atom-superposition-based method 



spin-polarized Tersoff-Hamann framework 



29 



3l] in the 



32] for simulating SP-STS by including the back- 



ground and tip-derivative terms into the calculated differential conductance following Pas- 
soni et al. [7(. We demonstrate that one can circumvent the numerical differentiation of the 
tunneling current for calculating dl/dV. The method is computationally cheap and it can be 
applied using results of any ab initio electronic structure code. Taking a prototype frustrated 
hexagonal antiferromagnetic system, a Cr monolayer on Ag(lll) in a noncollinear magnetic 
120° Neel state, we simulate differential conductance tunneling spectra and magnetic asym- 
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metries to illustrate the applicability of our method, and we analyze the contributing terms. 
The main advance of our tunneling model is the inclusion of the tip electronic structure, 



which is neglected in Refs. [32l . |33| , and it is only taken into account in a model way in Ref . 
Q| . Our method, based on first principles calculation of the tip electronic structure, enables 
to study tip effects on the SP-STS spectra. Note that a three-dimensional (3D) approach 
to STS has been presented recently, that is applicable to nonmagnetic systems only, and it 
takes into account the symmetry of the tip states but not the electronic structure of the 
tip apex 34J. Our model is also a 3D approach in the sense that we sum up contributions 
from individual transitions between the tip apex atom and each of the surface atoms assum- 
ing the one-dimensional (ID) Wentzel-Kramers-Brillouin (WKB) approximation for electron 
tunneling processes in all these transitions, thus we call it a 3D WKB approach. 

The paper is organized as follows: The theoretical model of SP-STS is presented in sec- 
tion [III As an application, we investigate the surface of one monolayer (ML) Cr on Ag(lll) 
in section IIII1 We simulate differential conductance tunneling spectra and magnetic asym- 
metries with two tip models, and we analyze the contributing terms to dl/dV. Moreover, we 
show an example for simulating two-dimensional (2D) differential conductance and magnetic 
asymmetry maps following a constant current contour above the surface, corresponding to 
a standard experimental setup. Our conclusions are found in section IIVI 

II. THEORETICAL MODEL OF SP-STS 

We start our model with the formulation of the tunneling current and the differential 
conductance in the framework of the one-dimensional (ID) WKB approximation, which 



will be used in our atom superposition approach later on. Assuming elastic tunne 
non-spinpolarized part of the tunneling current at zero temperature is given by 



ing, the 



E§,+eV 



I(V,d)=C T(E 1 V,d)n T (E)n s (E)dE J (1) 



E s 
F 



where V is the bias voltage, d the tip-sample distance, C an appropriate constant, Ep the 
Fermi energy of the sample surface, e the elementary charge, T the tunneling transmission 
coefficient, while riT{E) and ns{E) are the tip and sample densities of states, respectively. 
Performing a change of variable from E to U using E = Ep + ell, the tunneling current 
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reads 

t-V 



I(V,d) = Ce [ T(E s F + eU,V,d)n T {E s F + eU)n s {E s F + eU)dU. (2) 
Jo 

The applied bias voltage V in the tunnel junction defines the difference between tip and 
sample Fermi levels, E F = Ep + eV. Using this, the energy dependence of tit{E) can be 
rewritten related to the tip Fermi level Ep, and the tunneling current can be reformulated 

as 



r v 

I(V,d) = Ce T(E s F + eU,V,d)n T {El + eU -eV)n s (E s F + eU)dU. (3) 
Jo 

We denote the integrand by the formal quantity, 

^j(U, V, d) = CeT(E F + eU, V, d)n T (E T F + eU - eV)n s (E s F + eU), (4) 

called virtual differential conductance. The tunneling current can then be expressed as 

f v dl 

I(V,d) = —{U,V,d)dU. (5) 

The physical differential conductance can be obtained as the derivative of the tunneling 
current with respect to the bias voltage. This can formally be written as 



dl , Tr dl /rrrr f d dl / t j t ri i\ 



dU, (6) 

v=v 



— [T(Ep + eU, V, d)n T (E T F + eU - eV 



or using Eq.(j4]) as 

^(V, d) = CeT(E F + eV, V, d)n T {E T F )n s {Ep + eV) (7) 
«v 

n s (E b F + eU)dU. 

v=v 

This is a known formula in the literature [5|, [35|. If the tip electronic structure is assumed to 
be energetically flat, i.e. ut{E) = m, which is still a widely used approximation in the recent 
literature, then the ^/-dependence of riT(E F + eU — eV) disappears, i.e. riT(E F + eU — eV) = 
n F , and the differential conductance becomes 

^(V, d) = CeT(E F + eV, V, d)n T {E T F )n s {E s F + eV) (8) 
+ Cen T / — (£j + eU, V, d)n s (E s F + eU)dU 

JO 

Here, the second term is the so-called background term, which is a monotonous function of 



the bias voltage 



35l | . We would like to go beyond the assumption of the electronically flat 



tip by incorporating the tip electronic structure in the differential conductance expression. 
This way, the effect of the tip can be studied on the tunnelin g sp ectra. The explicit energy 
dependence of tit(E) can be calculated from first principles 22, 29|, or can be included in 
a model way jj. Following Eq.(J7J), the differential conductance can be reformulated as 

^(V, d) = CeT(E F + eV, V, d)n T {E T F )n s {E s F + eV) (9) 

+ Ce / — (Ef. + eU,V, d)n T {E T F + eU - eV)n s {E s F + eU)dU 
Jo 

+ Ce [ T(E S F + eU,V,d)^{E T F + eU - eV)n s (E s F + eU)dU. 
Jo 

Using that dn T (E F + eU — eV)/dV = —dn T (E F + eU — eV)/dU, the differential conductance 
at bias voltage V can be written as a sum of three terms, 

^(V, d) = ^j(V, V, d) + B(V, d) + D T (V, d) (10) 



with 

dl 



^j(V,V,d) = CeT(E s F + eV,V,d)n T (ET)n s (E s F + eV) (11) 

B(V,d) = Ce / —AE S F + eU,V,d)n T {E T F + eU - eV)n s {E s F + eU)dU (12) 
Jo vv 

D T {V,d) = -Ce J T{E s F + eU,V,d)^{E T F + eU - eV)n s (E s F + eU)dU. (13) 



Here, B(V, d) is the background term usually considered in recent STS theories By ,135), 
and D T (V,d) is a term containing the energy derivative of the tip density of states (DOS), 
which is rarely considered for practical STS calculations and analyses of experimental STS 
data. Including this term into our model, we can calculate the differential conductance 
naturally considering the effect of the tip electronic structure. Moreover, Eqs. ( fT0l -(fT3l) 
enable one to calculate STS spectra without the need for numerical differentiation of the 
tunneling current with respect to the bias voltage, thus, avoiding numerical errors. We 
demonstrate this in section II III for a complex magnetic surface probed with a magnetic tip. 
It can be shown that an alternative expression for the differential conductance can be derived 
using integration by parts, 

^(V, d) = ^(0, V, d) + B(V, d) + B 2 (V, d) - D S (V, d) (14) 
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with 



dl 

dU 



(0, V, d) = CeT(E F , V, d)n T {E T F - eV)n s {E s F ) 
dT 



(15) 



B 2 (V,d) = Ce / — {E S F + eU,V,d)n T {E T F + eU - eV)n s {E s F + eU)dU (16) 
Jo oU 

D s (V,d) = -Ce j T(E F + eU,V } d)n T (E T F + eU - eV)^-{E s F + eU)dU. (17) 

This way another background term, B2(V,d) enters the differential conductance formula, 
and the energy derivative of the sample DOS appears in the term D${V, d). The average of 
the two dl/dV expressions can also be formed as 



dI n/ ^ 1 
dV^ d)= 2 



§(o,v,d) + §(v,v,d) 



B{V, d) + l -B 2 (V, d) + ^ [D T (V, d) - D S (V, d)} , 

(18) 



which gives a third alternative form for the differential conductance within the ID WKB 
approximation. On the other hand, by subtracting Eq. (fT4|) from Eq. lflO]) . one gets 







[B 2 (V,d)-D T (V,d)-D s (V,d)]. 



(19) 



This is trivial since B2(V, d) — DT{V,d) — Ds{V,d) is related to the partial derivative of 
dI/dU(U, V, d) with respect to U: 



B 2 {V, d) - D T (V, d) - D S (V, d) = Ce 



v 



(20) 

From the three equivalent dl/dV formulas in Eqs. ffTOl) . ffT^j) and ( [18]) . the calculation of 
Eq.f riOl) needs the least mathematical operations, thus, we adopt this formula to our atom 



superposition approach 



29 



3l| in order to simulate STS spectra based on electronic structure 



data calculated from first principles. Note that the alternative expressions for dl/dV could 
prove to be useful for extracting information about the tip and sample electronic structures 
from experimentally measured spectra. 

Above derivations considered the non-spinpolarized part (TOPO) of the tunneling current 
and the differential conductance only. In magnetic STM junctions, the total tunneling 
current can be decomposed into non-spinpolarized (TOPO) and spin-polarized (MAGN) 



parts 



32 



33 



l36|,|37j, 



Itotal — Itopo + 1 



MAGN- 



(21) 
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Following the spin-polarized Tersoff-Hamann model 32[ and its adaptation to the atom 



superposition framework 



31 



331 ]. the magnetic contribution at a given energy is proportional 
to the scalar product of the tip and sample magnetic densities of states vectors, m T (E) and 
m s (E), respectively, 

^^L(E)^rn T (E)rn s (E). (22) 

Thus, the spin-polarized parts of dl/dV can similarly be calculated within the ID WKB 
approximation as reported above, just replacing nr(E)ns(E) by m T (E)m s (E) . 

Let us now formulate the tunneling current and the differential conductance and their 
TOPO and MAGN contributions within the atom superposition framework following Ref. 



311 ] . Here, we assume that electrons tunnel through one tip apex atom, and we sum up 
contributions from individual transitions between this apex atom and each of the surface 
atoms assuming the ID WKB approximation for electron tunneling processes in all these 
transitions. The tunneling current at the tip position R TIP (x, y, z) and at bias voltage V is 
given by 

Itotal(x, y, z, V) = I T opo(x, y, z, V) + I m agn{x, y, z, V), (23) 
where the contributions are given following Eq.([5]) as 

fV jr 

I T opo(x,y,z,V) = J -^°( x , y , z ,U,V)dU (24) 
Imagn(x, y, z, V) = J dlM ^ N (g» Vi z , U, V)dU. (25) 



Here, the integrands are the so-called virtual differential conductances 
Ttopo ( x v z U V) = e 2 - 



^^(x, y, z, U, V) = e 2 ^ ^ e - 2 « u ^ d "^n T {E T F + eU - eV)n a s (E s F + eU) (26) 



a 

-,2 



CLLJ It 

a 

respectively. Here, e 2 e 2 /h ensures that the dl/dU is correctly measured in the units of 
A/V. The sum over a corresponds to the atomic superposition and has to be carried out, 
in principle, over all surface atoms. Convergence tests, however, showed that including a 



relatively small number of atoms in the sum provides converged dl/dU values 29]. The tip 
and sample electronic structures are included into this model via projected DOS (PDOS) 
onto the atoms, i.e. ht{E) and ng(E) denote projected charge DOS onto the tip apex and 
the ath surface atom, respectively, while m T (E) and rn^{E) are projected magnetization 
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DOS vectors onto the corresponding atomic spheres. They can be obtained from collinear 



or noncollinear electronic structure calculations [31]. In the present work we determine the 
noncollinear PDOS for the sample surface and we use a collinear PDOS for a model CrFe 
tip 

Similarly to the tunneling current, the physical differential conductance can be decom- 
posed into non-spinpolarized (TOPO) and spin-polarized (MAGN) parts and it can be writ- 
ten at the tip position R TIP (x, y, z) and at bias voltage V as 

V, *, V) = ^(*, y, z, V) + ^^(z, y, z, V), (28) 

where the contributions are given following Eq. (|T0|) as 

**I2?°(x,v,z,V) = ^f^(x,y,z,V,V) + B TO Po(x,y,z,V) + D^ OPO (x,y,z,V) (29) 
^0^(x,y,z,V) = ^^(x,y,z,V,V)+B MAGN (x,y,z,V) + D^ AGN (x,y,z,Vm 

Obviously, the total differential conductance can also be written in the form of Eq.( llOj) as 

L (x, y, z, V) = ^-(x, y, z, V, V) + B(x, y, z, V) + D T (x, y, z, V) (31) 



dV K ' ' dU 



with 



§<x,y,z,V,V) = d J^fO( x , y , z ,V,V) + d -^( X ,y,z,V,V) (32) 
B(x,y,z,V) = B TO po(x,y,z,V) + B MAGN (x,y,z,V) (33) 
D T (x, y, z, V) = D™ PO (x, y, z, V) + D^ AGN (x, y, z, V). (34) 

The transmission probability for electrons tunneling between states of atom a on the surface 
and the tip apex is of the simple form, 

T{Ef + eU, V, d a (x, y, z)) = e -Mu,v)d a (x, y ,z) _ ^ 

This corresponds to a spherical exponential decay of the electron wavef unctions. Here 
d a (x, y, z) = \R TIP (x, y, z) — R a \ is the distance between the tip apex and the surface atom 
labeled by a with position vector R^. Assuming an effective rectangular potential barrier 
between the tip apex and each surface atom, the vacuum decay k can be written as 



«U,V) = U*n( ** + * + eV -eU), (36) 



where the electron's mass is m, h is the reduced Planck constant, and <ps and 0t are the 
average electron workfunction of the sample surface and the local electron workfunction of 
the tip apex, respectively. The method of determining the electron workfunctions is reported 



in Ref. 



311 ] . k is treated within the independent-orbital approximation |33l . 



3s 



which 



means that the same spherical decay is used for all type of orbitals. Extension of our model 



to take into account orbital dependent vacuum decay following Chen's work [40| is planned in 
the future, which is relevant for a more advanced description of tunneling between directional 
orbitals. Using Eq. (l35|) and the given form of the vacuum decay in Eq. (l36|) . the derivative 
of the transmission probability with respect to the bias voltage is obtained as 

dT fjps . tt rr A i \ \ me a ( \T(E$ + eU,V, d a (x,y, z)) 

gy{E F + eU,V,d a (x,y,z)) = -—d a (x,y,z) . (37) 

Considering these and following Eqs. (I12p and ffT3|) . the background and the tip-derivative 
contributions to dl/dV can be written as 

Trie 

B TOP o(x,y,z,V) = -e 2 -— (38 

ZTTfl 4 



rV -2 K (U,V)d a (x,y,z) 

V d a (x, y, z) / — — n T {E T F + eU - eV)n a s (E s F + eU)dU 

„ Jo k(U, V) 



x 

a 



me 3 



B MAGN (x,y,z,V) = -e — — (39) 

27in° 



rV -2K{U,V)d a (x,y,z) 

V d a (x, y, z) \ — — m T {E T F + eU - eV)m a s (E s F + eU)dU 

JO K \ U : V ) 



X 

a 



D™ PO (x,y,z,V) = -e* e - (40) 



2 



X 



J2 £ e - 2 < u ' v ^^^-{E T F + eU- eV)n%{E s F + eU)dU 



D^ AGN (x,y,z,V) = -e 2 j (41) 



x 



V- f V e -2 K (u,v) da (x,y,z)dmr^ E T + eU _ eV )rrf s (E s F + eU)dU. 
Jo 



Note that Eqs. (1161) and (IT7|) can similarly be formulated within the atom superposition 
approach, and two alternative expressions for the differential conductance can be given 
according to Eqs. (IT4T) and f|T5]) . 

I(x, y, z, V) and dI/dV(x,y,z,V) can be calculated at (x,y,z) grid points of a three- 
dimensional (3D) fine grid in a finite box above the surface. The recipe for simulating 
SP-STM images based on the 3D current map is given in Ref. j3l|. Here, we focus on the 
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simulation of dl/dV spectra. From the 3D differential conductance map, data can be ex- 
tracted that are directly comparable to experiments. For example, a single point spectrum 
corresponds to a fixed (xoiZ/oj^o) tip position, and two-dimensional (2D) spectra can also 
be obtained, where the image resolution is determined by the density of (x, y) grid points. 
There are usually two ways to define a 2D differential conductance map. The first method 
fixes the tip height at z = Zc = const and scans the surface, dI/dV(x,y, Zc,V), which is 
not frequently used in experiments since it does not ensure constant tunneling probability 
during the scan, and makes the interpretation of images difficult 15|. A constant tip-sample 
distance results in a constant tunneling transmission in the ID WKB model only. The sec- 
ond option measures dl/dV on a constant current contour, which is the widely used method 
in experiments. Simulation of this one can be done in two steps: First, we calculate the 
3D current map with the given bias voltage V s t a b, and at the second step we determine 
the height profile of a constant current contour, z(x,y,V sta b, Itotal — Ic — const), using 
logarithmic interpolation [3jJ. V sta b and I sta b = Ic are the tunneling parameters, and they 
stabilize the tip position at the height of z(x, y, V sta b, Itotal = Istab) above the (x, y) sam- 
ple surface point. The differential conductance map on the constant current contour is then 
given by dl /dV(x,y, z(x,y,V sta b, Itotal = Istab), V), where the ^/-dependence is obtained 
by sweeping the bias voltage range using a lock-in technique in experiments However, 
since the total current also contains a magnetic contribution, this method does not provide 
dl/dV measurement at a constant tunneling probability due to the enhanced corrugation 
of the magnetic image. Spin-resolved spectroscopic mode of SP-STM [jj]] can correctly be 
simulated by replacing Itotal in the previous formula with Itopo, that truly ensures con- 
stant tunneling transmission while the spectrum dI/dV(x, y, z(x, y, V sta b, Itopo — Istab), V) 
is taken. This quantity can easily be calculated within the framework of our 3D WKB 
method, while in experiments this is more difficult to achieve since one can directly control 
the total tunneling current during the measurement, which contains both Itopo and Imagn 
contributions [is! ]. 

By simulating differential conductance spectra above a magnetic surface with parallel 
(P) and antiparallel (AP) tip magnetization directions with respect to the magnetization 
direction of a given magnetic surface atom, the so-called magnetic asymmetry can be defined 



2l| . In our case this quantity can be calculated at all considered positions of the tip apex 
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atom, i.e. at all (x, y, z) grid points within our finite box above the surface: 

dI p /dV(x, y z V) - dI AP /dV(x, y z V) 
A MAGN (x, y, z, V) = dlP/dvix ,y, z ,v) + dI">/dV{x,y,z,vy ^ 

From this, the magnetic asymmetry can similarly be calculated on constant current contours 
as described in the previous paragraph. We will demonstrate this above a complex magnetic 
surface in the next section. Using Eq.( l28l) . and the fact that the magnetic contribution for 
the AP tip magnetization direction dI M p AGN / dV equals —dI MAGN /dV, since the tip magne- 
tization DOS vector m T (E) changes sign at all energies compared to the P tip magnetization 
direction, the differential conductances take the following form: 

dl p / dV(x, y, z, V) = dI TO po/ dV(x, y, z, V) + dl[ IAGN /dV(x, y, z, V) 
dI AP /dV(x, y, z, V) = dI TO Po/dV(x, y, z, V) - dI p AGN /dV(x, y, z, V). (43) 

Thus, the magnetic asymmetry can be expressed as the fraction of the MAGN and TOPO 
differential conductances from Eqs. (|2"9l) and (1301) as 



A dI/dV (r v z V) - dI MAGN/dV(x,y,z,V) 

A MAGN (x, y, z, V) - dI P opo/dV{Xj ^ z v) (44) 

mLagn/ du (g» y, g> Yi v ) + b magn(. x , y, g> Y) + d t iagn ' p ( x > y» z , Yl 

dI TO Po/dU(x, y, z, V, V) + B TO po(x, y, z, V) + D™ p °(x, y, z, V) 

This is the correct magnetic asymmetry expression based on the physical differential con- 
ductances. However, a magnetic asymmetry can similarly be defined taking the virtual 
differential conductances from Eqs. ( 1261) and ( 1271) : 

A di/du ( rr\ dI p AGN /dU{x,y,z,V,V) 

/i MAGN\ x i Vi Z > v ) ~ 777 73777 t 7 77\ ■ l 40 J 

dlroPo/dU [x, y, z, V, V) 

There are even more possibilities to define magnetic asymmetries, by adding the background 
terms in Eqs. (138j) and ( 1391) . or the tip-derivative terms in Eqs. ( HU1) and (14 lj) to the corre- 
sponding virtual differential conductance and then performing the division: 

A dI/dU+B, , trt dI MAGN/ dU fa Vj *, Y V ) + B MAGN( X > Vj z > V ) ( a ^ 

magn & V,z,V)- dlTOpo/dU{Xj y> Zj Vi v) + Btopo {x, y, z, V) ' l4bj 



A di/du+D T , dI^ AGN /dU (x, y, z, V, V) + Dt AGN,p (x, y, z, V) 

magn [x,y,z,V)- dlTOpo/du{x ^ iZiViV) + D TOPO {x ^^ v) • (47) 



As dI/dU(V, V) is one component of dI/dV(V) according to Eq. (l3Tj) . we will compare them 
and also the magnetic asymmetry expressions in Eqs. (l44p - (l4"Tj) . in order to estimate the 
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error one makes when neglecting the background and tip- related components of dI/dV(V) 
for a given combination of a complex magnetic surface and a magnetic tip. 

It has to be noted that the presented method can also be applied to study nonmagnetic 
systems, where all magnetic contributions equal zero and the corresponding topographic 
STS spectra can be simulated. Of course, in this case, the magnetic asymmetry is zero. 



III. RESULTS AND DISCUSSION 



In order to demonstrate the reliability and capabilities of our model for simulating SP- 
STS on complex magnetic surfaces, we consider a sample surface with noncollinear magnetic 
order. One ML Cr on Ag(lll) is a prototype of frustrated hexagonal antiferro magnets 33]. 
Due to the geometrical frustration of the antiferromagnetic exchange interactions between 
Cr spin moments, its magnetic ground state is a noncollinear 120° Neel state 32]. In the 
presence of spin-orbit coupling, two types of Neel states with opposite chiralities can form, 
and one of them is energetically favored 31]. 

We performed geometry relaxation and electronic structure calculations based on Density 
Functional Theory (DFT) within the Generalized Gradient Approximation (GGA) imple- 
mented in the Vienna Ab-initio Simulation Package (VASP) 
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431 ] . A plane wave basis set 



for electronic wavefunction expansion together with the projector augmented wave (PAW) 



method 



44j has been applied, while the exchange-correlation functional is parametrized ac- 



cording to Perdew and Wang (PW91) 45|. For calculating the fully noncollinear electronic 



structure we used the VASP code as well 
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471 ]. with spin-orbit coupling considered. 



We model the Cr/Ag(lll) system by a slab of a five-layer Ag substrate and one monolayer 
Cr film on each side, where the surface Cr layers and the first Ag layers underneath have 
been fully relaxed. A separating vacuum region of 14.6 A width in the surface normal 
(z) direction has been set up between neighboring supercell slabs. The average electron 
workfunction above the surface is <ps — 4.47 eV. We used an llxllxl Monkhorst-Pack 
(MP) [48( k-point grid for calculating the pro ject ed electron DOS onto the surface Cr atoms 
in our (a/3 x a/3) magnetic surface unit cell 31]. 



Dependence of the tunneling spectra on the tip electronic structure can be studied by 
considering different tip models. In this work we compare spectra and magnetic asymmetries 
measured by a magnetic CrFe tip and an electronically flat magnetic tip. The electronic 
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structure data of the CrFe tip was taken from Ref. [22], where the tip was modeled as a 
single Cr apex atom on the Fe(OOl) surface. Ferriani et al. furthermore reported that an 
antiferromagnetic coupling of the Cr adatom to the Fe(OOl) surface is energetically preferred, 
and the vacuum spin polarization is fairly constant at around +0.8 in the energy range 
[Ep — leV, Ep + leV] 22j. The local electron workfunction above the tip apex is assumed to 
be 4>t = 4.5 eV, that has been used to obtain the energy dependent vacuum decay in Eq. (|36|) . 
The electronically flat magnetic tip has been modeled based on the electronic structure of 
the Cr apex of the CrFe tip. The charge and absolute magnetization DOS, tit(E) and 
\i71t{E)\, respectively, have been averag ed in the [E% - 2eV, E T F + 2eV] range. We obtained 
rip = 1.33/ eV and nip = 1.06/eV, and, thus, the spin polarization is Pt = rap/nT = +0.8. 
In this case, the tip-derivative term of the differential conductance Dt(V) is zero, since 
driT{E)/dE = drriT{E)/dE = 0. The vacuum decay can be modeled using Eq. (l36j) . where 
k(U, V) has an explicit ^/-dependence, and we assume that 4>t — <Ps- Alternatively, a simpler 
model for k(U) can be considered without ^-dependence: 

k(U) = \^2m (<P s -eU). (48) 

In this case the background term of the differential conductance B(V) is zero, since the 
tunneling transmission does not depend on the bias voltage, and the physical differential 
conductance equals the virtual differential conductance, i.e. dI/dV(V) = dI/dU(V,V). On 
the other hand, by assuming a ^-dependent vacuum decay k(U,V), B(V) is not zero and 
it contributes to the total differential conductance, i.e. dI/dV(V) = dI/dU(V, V) + B(V). 

Figure [1] shows the bias dependence of the total tunneling current I{V), calculated using 
Eq. (l23i) . at the position z = 3.5 A above a surface Cr atom probed with the CrFe tip having 
parallel (P) magnetization direction compared to the underlying surface Cr atom. Positive 
current means tunneling from the tip to the sample surface, whereas the current is nega- 
tive in the opposite direction. We find that the absolute value of the current is higher in 
the negative bias range compared to the positive range. This is due to the surface and tip 
electronic structures. The sample occupied DOS combined with the tip unoccupied DOS is 
greater than the sample unoccupied DOS combined with the tip occupied DOS. Performing 
a numerical differentiation of /(V) with respect to V, we obtain the differential conductance 
at this particular tip position. As can be seen this is extremely noisy, and a smoothing 
procedure should be applied to it before further analysis. Alternatively, the differential 



14 



conductance can be calculated using Eq.( l3Tj) . implemented within the atom superposition 
approach. Figure [1] shows that dl/dV obtained this way (black curve) is a smooth function 
that fits precisely to the noisy numerical derivative of the current. Thus, we demonstrated 
that one can circumvent the numerical differentiation of the tunneling current and a sub- 
sequent smoothing for calculating dl/dV based on first principles electronic structure data. 
We obtain more information about the dl/dV by analyzing its components, the virtual dif- 
ferential conductance dI/dU(V, V), the background term B(V), and the tip- derivative term 
Dt(V). One can recognize that most part of the dJldV peak structure is already included 
in the dI/dU(V,V) term. As discussed in Ref. 29|, the essential features of the virtual 
differential conductance are formed by matching the sample and tip electronic structures, 
i.e. the bias voltage fixes the relative energetic position of tip and sample electron states. 
Apart from this, the peak structure of Dt{V), calculated via Eqs. P0|) - fj4T|) . clearly shows 
up in the dl/dV, particularly pronounced at high bias voltages. The reason is the rapidly 
changing tip electronic structure in these energy regions. The features from dI/dU(V, V) 
and Dt(V) are transferred to the physical differential conductance, since B(V), calculated 
via Eqs. ( I38i) -( l39l) . is smooth compared to the other two components in the whole bias range. 
Moreover, we find that B(V) is a monotonous function of the bias voltage, and it is ne arly 
proportional to I(V) as has been reported earlier for different levels of STS theories Q] • 
The proportionality function B(V)/ I(V) is plotted in the inset of Figure [TJ It can be seen 
that its sign is in agreement with Ref. [3] and it has a non-trivial bias dependence. The 
B(V)/I(V) function could, in principle, be calculated at different tip-sample distances (z), 
and could be compared to analytical expressions denoted by f(z, V) reported in The 
comparison is, however, not straightforward due to two reasons. First, the analytical ex- 
pressions were reported based on the ID WKB approximation, whereas our model is a 3D 
atomic superposition approach based on WKB, which results in an effective transmission 
coefficient, different from the ID WKB transmission. Note that a 3D approach to STS with 
another effective transmission coefficient has recently been reported by Donati et al. 34]. 
Second, in Figure Q] we reported the sum of the TOPO and MAGN contributions, while 
the related STS literature is concerned with nonmagnetic systems only, which corresponds 
to the analysis of the topographic part of the spin-polarized results. Consideration of the 
spin-polarized tunneling complicates the analytical calculations that are unavailable at the 
moment. The analysis of B(z, V) / I(z, V) along the discussed lines could be a future research 
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direction that is beyond the scope of the present study. In the following we focus on the 
comparison of SP-STS spectra by reversing the tip magnetization direction, and also using 
the flat magnetic tip model. 

Figure [2] shows simulated single point differential conductance tunneling spectra following 
Eq.f l3ip . probed with the flat magnetic tip and the model CrFe tip, z = 3.5 A above a surface 
Cr atom. Parallel (P) and antiparallel (AP) tip magnetization directions are set relative 
to the underneath surface Cr atom. It can clearly be seen that measuring the spectra 
with oppositely magnetized tips of the same type result in different differential conductance 
curves, in agreement with SP-STS experiments performed on oppositely magnetized sample 



with a fixed tip magnetization direction 
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2l|. For the flat magnetic tip, two different 



vacuum decays, n{U) and k(U, V) are assumed using Eqs. (jUJ) and (13T)|) . respectively. For 



ow 



29| 



the bias-independent vacuum decay (dotted curves) we find that dI p /dV > dI AP /dV be 
V = +0.54 V, while dI p /dV < dI AP /dV above V = +0.54 V. In our previous work 
we identified the effective spin polarization [P T (E)P S (E) = m T (E)m s (E) / (riT(E)ns(E))} 
responsible for this effect. This is the decisive factor for determining the sign of the magnetic 
contribution to dl/dV at energy E in the improved SP-STS model presented in section HT1 
as well. The magnetic part of the physical differential conductance is given in Eq.( |30|) . 
Since the vacuum decay does not depend on the bias voltage V for the dotted curves, and 
the tip is electronically flat, dluACN / dV{V) = dlMAGN I 'dU(V,V). Thus, the sign change 
of dlMAGN I dV occurs at the sign change of dlMAGN / dU{V,V), i.e. at the reversal of the 
sample spin polarization vector at 0.54 eV above the sample Fermi level 3l|. For the flat 
magnetic tip and the assumed bias dependent vacuum decay (dashed curves) we find that 
dI p /dV > dI AP /dV below V = +0.5 V, and dI p /dV < dI AP /dV above V = +0.5 V, 
i.e. the sign change of the magnetic component is slightly shifted toward zero bias. The 
reason is the nonzero background term Bmagn(V) due to k(U,V), and dlMAGN / dV(V) = 
dI MAGN /dU(V,V) + B MAGN (V) has to be considered. Note that D^ AGN (V) is still zero 
because of the constant tip magnetization DOS. Comparing the two vacuum decay models 
for the flat tip, it is clear that the topographic part of the background term has another 
effect on the heights of the spectra, i.e. they are enhanced and reduced in the negative and 
positive bias ranges, respectively, compared to the k(U) model. On the other hand, the 
features of the spectra (peaks and dips) occur at the same bias positions for both vacuum 
decay models. 



1(3 



The inclusion of a realistic tip electronic structure into our model complicates the spec- 
tra even more. This is demonstrated in Figure [2] for the CrFe tip model (solid lines). 
In this case all three terms contribute to the differential conductance, and dI/dV(V) = 
dI/dU(V,V) + B(V) + D T (V). Thus, the relative heights of the differential conductance 
tunneling spectra dI p /dV and dI AP /dV are determined by the superposition of the mag- 
netic dluAGN / dU{V, V), Bmagn(V), and D^ IAGN (V) terms. The role of the effective spin 
polarization is more complicated, since, apart from the dlMAGN / 'dU(V,V) term, it appears 
in the dl/dV expression through the bias-integrated quantities Bmagn(V) and D^ IAGN (V). 
For the P tip magnetization, dI p /dV is the same as the black solid curve in Figure dj and 
its contributions are also shown there. In Figure [21 we observe more changes of the relative 
height of the dI p /dV and dl AP jdV spectra measured with the CrFe tip than with the flat 
tip. These include the sign changes of the magnetic part of the spectra, similarly as before. 
We find that dI p /dV > dI AP /dV in the bias interval [-1.04 V, +0.49 V], and a reversed 
relation is obtained in the complementary bias regime. Comparing the spectra to the ones 
measured with the flat magnetic tip, we see that they are qualitatively closer to the k(U, V) 
model used for the flat tip due to the presence of the background terms. Moreover, the indi- 
vidual features coming from the sample and the tip electronic structures can be identified. 
It would be quite helpful analyzing experimental tunneling spectra following the recipe: 

(1) Measure dl/dV spectra experimentally. 

(2) Simulate dl/dV spectra on a reliable model surface with an electronically flat tip. 

(3) Simulate dl/dV spectra on the sample surface with a reliable model tip corresponding 
to the experiment. 

(4) Identify the tip- and sample-related peaks. 

This procedure could be applied for both nonmagnetic and magnetic STS setups. In the 
magnetic case, the spin polarization of the flat tip can be tuned to change the magnetic 
sensitivity and identify surface magnetic contributions. In our case we identify the peak at 
-1.2 V, indicated by a vertical dotted line in Figure EJ coming from the CrFe tip electronic 
structure since it is missing from the spectra calculated with the flat tip. 

The relative heights of the differential conductance tunneling spectra dI p /dV and 
dl AP jdV can also be determined from the magnetic asymmetry, Eq. (f42|) . Let us compare the 
magnetic asymmetries calculated from the spectra in Figure |2] using the two magnetic tips. 
Moreover, for the CrFe tip we compare the Amagn expressions defined in Eqs. (l4~4"l) - (14T1) . in 



17 



order to estimate the error one makes when neglecting the background and tip-related com- 
ponents of dI/dV(V). Figure [3] shows the calculated Amagn(V) functions at z — 3.5A above 
a surface Cr atom. It can be seen that A^agjv (V) and A^ac/n (V) (dashed curves with 
triangle symbol) behave qualitatively similarly. In addition, ^2^(7) is greater than 
^magn (^0 m almost the full studied bias range. The opposite relation holds between 
V and +0.3 V only, however, the relative difference between the two quantities is less than 
1.4 % in this regime. Moreover, these two magnetic asymmetries are within 5% relative 
difference in the bias range [-0.23 V, +0.31 V]. 

Considering the CrFe tip, ^^'^^(^l (black solid curve) is qualitatively different from 
the two A MAGN (V) functions calculated with the flat tip, e.g. it has a richer structure at 
positive bias voltages. More importantly, there is an extra sign change of A^Jg^ r//dV (V) 
occurring at -1.04 V apart from +0.49 V. These correspond to the height changes of dl p jdV 
and dl AP I dV relative to each other in Figure [2j Let us estimate the error of Amagn when ne- 
glecting the background and the tip-derivative terms. According to Eq.( l45"I) . A^^q^ L (V) 
(curve with symbol 'O') considers the virtual differential conductances only. It is within 10% 
relative error compared to A < l^Q^ dV (V) in the bias range [-0.65 V, +0.1 V]. However, its 
sign does not correspond to A^^q^^ (V) in the bias intervals [-2 V, -1.04 V] and [+0.49 
V, +0.54 V]. Adding the background term B{V) to dI/dU(V, V) results in an improved 
differential conductance expression, and A^^Q^^ dU+B (V) (curve with symbol '+'), defined 
in Eq. pS"]) . behaves qualitatively similarly to A c ^^^ dh \V) above -0.65 V. However, its sign 
change is shifted to +0.45 V from +0.54 V. On the other hand, a sign change in the negative 
bias range appears similarly to A c h ^Q^ dV (y), though it is shifted to -1.62 V. In the low 
bias range, A ( ^^Q^ dU+B (V) is within 10% relative error compared to A^^gn ^ '(V) * n a 
decreased bias range of [-0.34 V, +0.1 V]. Finally, by adding the tip-derivative term Dt(V) 
to dI/dU(V, V), A^Jagn T ^ dU+Dl \V) (curve with symbol 'x'), defined in Eq. pT)) . shows the 
most closely related shape to A^^gn V (V)- Furthermore, it is also quantitatively close to 
the physical magnetic asymmetry as its sign changes occur at -1.01 V and +0.5 V, and it is 
within 10% relative error compared to A c ^^^ dV (V) in an increased bias interval [-0.90 V, 
+0.45 V]. Summarizing this paragraph, the contribution of all three terms to the dI/dV(V) 
according to Eq. (1311) is needed to define the physical magnetic asymmetry, that should be 
comparable to experiments. 

Our method presented in section [TT] also enables one to simulate two-dimensional (2D) 
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dl/dV and magnetic asymmetry maps in high spatial resolution above the surface, that can 
be compared to results of SP-STS experiments. Such experiments are supposedly performed 
at a constant tunneling probability 15], i.e. the tip scans the surface following a constant 
TOPO current contour. Figure H] illustrates this capability of our method, where we used 
the flat tip model. By choosing V sta b = +1 V, we calculate the 3D TOPO current map in the 
box above the surface shown in the top left part of Figure HI From this 3D data we extract 
the current contour of Itopo — 0.055 nA, which is around 3.5 A above the sample surface 
and has a corrugation of 3.5 pm. This contour is plotted inside the box in the top left part 
of Figure HI The surface scan area and the magnetic unit cell are shown in the top right part 
of Figure HI indicated by a black-bordered rectangle, and a yellow (light gray) rhombus, 
respectively. For calculating the 2D dl/dV map, we vary the z position of the tip apex 
atom following the current contour above each (x, y) lateral points. By considering P and 
AP tip magnetization directions relative to the (1/2, \/3/2) direction (i.e. the magnetization 
direction of the surface Cr atom at the bottom left corner of the scan area, see top right 
part of Figure HI), spin-resolved dl/dV and magnetic asymmetry maps can be obtained at 
different V bias voltages. As an example, dl p /dV(x,y) and A d ^Q N (x,y) are shown in the 
bottom part of Figure H] calculated at V = +1 V. At this bias, all surface Cr spin polarization 
vectors point opposite to their local magnetic moment directions 3l|, and since Pt = +0.8 
is set with respect to the (1/2, y/3/2) direction, the leading term of the magnetic differential 
conductance, dluAGN / dUiV, V) is negative above the surface Cr atom at the bottom left 
corner of the scan area. Moreover, the sign of dlyiAG n / dU (V, V) changes to positive above 
the other two Cr atoms in the magnetic unit cell. These, together with the fact that the 
topographic differential conductance is the same above all Cr atoms, result in the minimal 
total dl p /dV(x,y) above the Cr atom at the bottom left corner of the scan area (0.074 
nA/V, magn. moment parallel to tip), whereas above the other two Cr atoms dl/dV is 
maximal (0.102 nA/V, magn. moment not in line with tip). The indicated constant dl p jdV 
contours illustrate this finding as well. For the same reason, the AP tip magnetization 
direction would result in a reversed dI AP /dV map concerning the heights above the non- 
equivalent magnetic Cr atoms, i.e. it would be maximal above the Cr atom at the bottom 
left corner of the scan area, and minimal above the other two Cr atoms. Thus, the difference 
between dl p /dV(x, y) and dI AP /dV(x, y) is minimal and negative above the bottom left Cr 
atom in the scan area, and maximal and positive above the other two Cr atoms. These 
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findings explain qualitatively well the simulated Amagn(x, y) map, where Amagn has a 
minimum above the Cr atom at the bottom left corner, and maximum above the other 
;wo Cr atoms. Since the effective spin polarization (ESP) is proportional to A MAGN (x,y) 



2l| . a qualitatively correct ESP(x, y) map above the surface has been simulated within the 
indicated approximations. The Amagn(%, y) — contour corresponds to zero ESP, which 
acts as a border between surface regions with positive and negative ESP at the given bias. 
Note that the sign of the tip spin polarization has a crucial effect on the ESP(x, y) map. 
Reversing the sign of Pt would result in a reversed ESP(x, y) map. 

We suggest that by applying our method to magnetic surfaces, two-dimensional 
dl p /dV(x,y), dI AP /dV(x,y), A MAGN (x,y), and, thus, a qualitative ESP(x,y) map can 
be constructed on arbitrary TOPO current contours (i.e. at different tip-sample distances), 
and at arbitrary V bias. By explicitly considering the tip electronic structure in our SP-STS 
model based on experimental information, it would help in a more reasonable interpretation 
of experimentally measured tunneling spectra and magnetic asymmetries. 

IV. CONCLUSIONS 



We extended our atom-superposition-based method [3l| for simulating spin-polarized 
scanning tunneling spectroscopy by including the background and tip-derivative terms into 
the differential conductance formula following Passoni et al. We demonstrated that one 
can circumvent the numerical differentiation of the tunneling current for calculating the 
dl/dV. Furthermore, we showed that our simulated data can be related to standard experi- 
mental setups. Taking the tip electronic structure into account, the effect of a richer variety 
of electronic structure properties can be investigated on the tunneling transport within the 
indicated approximations (atom superposition, spherical vacuum decay). The method is 
computationally cheap and it can be applied based on the results of any ab initio electronic 
structure code. Taking a prototype frustrated hexagonal antiferromagnetic system, a Cr 
monolayer on Ag(lll) in a noncollinear magnetic 120° Neel state, we simulated differential 
conductance tunneling spectra and magnetic asymmetries to illustrate the applicability of 
our method, and we analyzed the contributing terms. We found that the features of the 
tunneling spectra are coming from the virtual differential conductance and tip-derivative 
terms, and the background term is proportional to the tunneling current. We showed evi- 
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dence that the tunneling spectra and the related magnetic asymmetries are sensitive to the 
tip electronic structure and to the vacuum decay. We also simulated a 2D dl/dV map and 
a qualitatively correct effective spin polarization map of the magnetic surface following a 
constant current contour. In order to help extracting information about the tip and sam- 
ple electronic structures from experimentally measured tunneling spectra, we derived three 
equivalent expressions for the dl/dV within the ID WKB approximation, and we suggested 
a combined experimental-theoretical procedure to identify tip- and sample-related features. 
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* *B(V)/I(V) 




Bias Voltage [V] 



FIG. 1: (Color online) Comparison of single point differential conductance tunneling spectra calcu- 
lated from numerical differentiation of the tunneling current /(V), and dl/dV calculated according 
to Eq. (|3ip . and its contributing terms, the virtual differential conductance dI/dU(V, V), the back- 
ground term B(V), and the tip-derivative term Dt(V). The model CrFe tip apex is 3.5 A above 
a surface Cr atom and its magnetization direction is parallel to that of the underlying surface Cr 
atom. The inset shows the ratio of B(V)/I(V). 
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FIG. 2: (Color online) Comparison of simulated single point differential conductance tunneling 
spectra following Eq. (|3ip . probed with the flat magnetic tip and the model CrFe tip, 3.5 A above 
a surface Cr atom. Parallel (P) and antiparallel (AP) tip magnetization directions are set relative 
to the underneath surface Cr atom. For the flat magnetic tip, two different vacuum decays, n(U) 
and k(U,V) are assumed using Eqs. (JUJ) and ([36]) . respectively. The vertical dotted line at -1.2 
V shows the bias position of the identified STS peak coming from the electronic structure of the 
CrFe tip. 
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FIG. 3: (Color online) Comparison of magnetic asymmetries 3.5 A above a surface Cr atom probed 
with the flat magnetic tip and the model CrFe tip. A^ a ^Q^\ A^^^^'^, and ^magtv^^ are 
calculated from the corresponding P and AP spectra shown in Figure [2J For the CrFe tip we 
compare the Amagn expressions defined in Eqs. (|i^) - (|17|) . 
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FIG. 4: (Color online) Top left: The three-dimensional finite box for our data collection above 
the Cr/Ag(lll) surface, containing a constant current contour about 3.5 A above the surface with 
iTOPoiy = +1V)=0.055 nA calculated with the flat magnetic tip. Top right: Surface geometry 
of 1 ML Cr on Ag(lll). The Cr and Ag atoms are denoted by spheres colored by green (medium 
gray) and purple (dark gray), respectively, while the magnetic moments of individual Cr atoms are 
indicated by (red) arrows. The (V3 x y/3) magnetic unit cell is drawn by yellow (light gray) color. 
Bottom: Simulated 2D differential conductance map, dl p /dV(x,y, z(x,y,V = +1V, Itopo = 
0.055nA),y = +1V), and magnetic asymmetry map, A d ^/^ N (x,y, z(x,y,V = +1V, Itopo = 
0.055nA),y = +1V), using the flat tip. Black lines correspond to constant value contours. The 
color coding is given in the top left part of the figure. The surface Cr positions are denoted by 'X'. 
During the scan the vertical tip position z(x, y) follows the constant current contour shown in the 
top left part of the figure. Parallel (P) and antiparallel (AP) tip magnetization directions are set 
relative to the surface Cr atom at the bottom left corner of the scan area. 
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